close all
clear
warning off

%% scenario
load('location_1202.mat');
fprintf('%d cells and %d users per cell\n',N,K);
plot_loca(loc_BS_ATG,loc_airplane,loc_BS,loc_user,num_terres_BS_layer,dist_BS_half);

MA=8;
MG=8;

%% initialize noise
% noise_bandwidth = 180*1000; % 180kHz
% noise_bandwidth = 10*10^6; % 10MHz
noise_bandwidth = 0.5*10^6; % 0.5MHz
noise_temp = 300;
Noise_power_spectral_densityWHz = 1.38064852*10^(-23)*noise_temp; % W/Hz
Noise_power_spectral_densitymWHz = 1000*Noise_power_spectral_densityWHz; % mW/Hz
Noise_power_spectral_densitydBmHz = 10*log10(Noise_power_spectral_densitymWHz); % dBm/Hz
% Noise_power_spectral_densitydBmHz = -170; % dBm/Hz
noise = 10^(Noise_power_spectral_densitydBmHz/10)*noise_bandwidth; % mW!!!!!!!!!!!!!!!!!!!!!!!!!!!
% noisedBm = 10*log10(noise); % =-117.4473dBm
% noisedBm = -170+10*log10(180*1e3); % =-117.4473dBm
% 对噪声归一化 %???????????????????????????
% Weighted Sum-Rate Maximization for Reconfigurable Intelligent Surface Aided Wireless Networks


%% constraints ???????????????????????????????????????????
% maximum interference power of satellite terminal (mW) 
% interfere_max=2*10^(-9);
interfere_max=1*10^(-8);
% interfere_max=1*10^(-14);

% maximum transmit power of TBS(mW)
pow_max=60*1000*ones(N,1);
% maximum transmit power of ABS(mW)
% pow_air=60*1000;
% xlabel_set = [80,90,100];


%% parameters
paraSystem.MA=MA;
paraSystem.MG=MG;
paraSystem.N=N;
paraSystem.K=K;
paraSystem.freq_=4*10^9; % 4GHz C band
% paraSystem.freq_=18*10^9; % 18GHz Ka band
paraSystem.noise_bandwidth=noise_bandwidth;
paraSystem.Noise_power_spectral_densitydBmHz=Noise_power_spectral_densitydBmHz;
% paraSystem.interfere_max=interfere_max;
paraSystem.pow_max=pow_max;
% paraSystem.pow_air=pow_air;
% paraSystem.SINR_min=SINR_min;



%% Generate the channel coefficients
[hGA,hGG] = generate_channel1127(paraSystem,loc_airplane,loc_BS,loc_user);
% 计算信道sqrt(h'*h)的数量级
[avechannel,count_inter,count_intra]=channel_analysis(N,K,hGA,hGG);
%         hhhhhhhhhhhhh

%% Direct FP
[ww_D_FP,auxi_vari,num_loop] = alg_Direct_FP(N,K,MG,hGG,noise,pow_max);
cons_sum_power_D_FP = cal_sum_power(N,K,ww_D_FP);
fprintf('%d/%d: the transmit power constraints of each terrestrial BS\n',...
    sum((pow_max- cons_sum_power_D_FP)>-sum(pow_max)/length(pow_max)*1e-3),length(pow_max));
cons_sum_interfe_D_FP = cal_sum_interfe(N,K,ww_D_FP,hGA);
fprintf('%d/1: the interference power constraint of the primary network\n',...
    (interfere_max- cons_sum_interfe_D_FP)>-interfere_max*1e-3);
sum_rate_terres_D_FP = cal_sum_rate_terres(N,K,ww_D_FP,hGG,noise);
fprintf('sum_rate_terres (alg_Direct_FP) = %d (bps/Hz)\n',sum_rate_terres_D_FP);

%% Direct FP in cognitive network
[ww_D_FP_cog,auxi_vari_cog,num_loop_cog] = alg_Direct_FP_cog(N,K,MG,hGG,hGA,noise,pow_max,interfere_max);
cons_sum_power_D_FP_cog = cal_sum_power(N,K,ww_D_FP_cog);
fprintf('%d/%d: the transmit power constraints of each terrestrial BS\n',...
    sum((pow_max- cons_sum_power_D_FP_cog)>-sum(pow_max)/length(pow_max)*1e-3),length(pow_max));
cons_sum_interfe_D_FP_cog = cal_sum_interfe(N,K,ww_D_FP_cog,hGA);
fprintf('%d/1: the interference power constraint of the primary network\n',...
    (interfere_max- cons_sum_interfe_D_FP_cog)>-interfere_max*1e-3);
sum_rate_terres_D_FP_cog = cal_sum_rate_terres(N,K,ww_D_FP_cog,hGG,noise);
fprintf('sum_rate_terres (alg_Direct_FP in cognitive network) = %d (bps/Hz)\n',sum_rate_terres_D_FP_cog);


%% WMMSE
[ww_WMMSE,uu,vv,num_loop] = alg_WMMSE(N,K,MG,hGG,noise,pow_max);
cons_sum_power_WMMSE = cal_sum_power(N,K,ww_WMMSE);
fprintf('%d/%d: the transmit power constraints of each terrestrial BS\n',...
    sum((pow_max- cons_sum_power_WMMSE)>-sum(pow_max)/length(pow_max)*1e-3),length(pow_max));
cons_sum_interfe_WMMSE = cal_sum_interfe(N,K,ww_WMMSE,hGA);
fprintf('%d/1: the interference power constraint of the primary network\n',...
    (interfere_max- cons_sum_interfe_WMMSE)>-interfere_max*1e-3);
sum_rate_terres_WMMSE = cal_sum_rate_terres(N,K,ww_WMMSE,hGG,noise);
fprintf('sum_rate_terres (WMMSE) = %d (bps/Hz)\n',sum_rate_terres_WMMSE);

%% Closed Form FP
[ww_CF_FP,uu,vv,num_loop] = alg_Closed_Form_FP(N,K,MG,hGG,noise,pow_max);
cons_sum_power_CF_FP = cal_sum_power(N,K,ww_CF_FP);
fprintf('%d/%d: the transmit power constraints of each terrestrial BS\n',...
    sum((pow_max- cons_sum_power_CF_FP)>-sum(pow_max)/length(pow_max)*1e-3),length(pow_max));
cons_sum_interfe_CF_FP = cal_sum_interfe(N,K,ww_CF_FP,hGA);
fprintf('%d/1: the interference power constraint of the primary network\n',...
    (interfere_max- cons_sum_interfe_CF_FP)>-interfere_max*1e-3);
sum_rate_terres_CF_FP = cal_sum_rate_terres(N,K,ww_CF_FP,hGG,noise);
fprintf('sum_rate_terres (alg_Closed_Form_FP) = %d (bps/Hz)\n',sum_rate_terres_CF_FP);


% %% Closed Form FP in cognitive network
% [ww_CF_FP_cog,uu,vv,num_loop] = alg_Closed_Form_FP_cog1(N,K,MG,hGG,hGA,noise,pow_max,interfere_max);
% cons_sum_power_CF_FP_cog = cal_sum_power(N,K,ww_CF_FP_cog);
% fprintf('%d/%d: the transmit power constraints of each terrestrial BS\n',...
%     sum((pow_max- cons_sum_power_CF_FP_cog)>-1),length(pow_max));
% cons_sum_interfe_CF_FP_cog = cal_sum_interfe(N,K,ww_CF_FP_cog,hGA);
% fprintf('%d/1: the interference power constraint of the primary network\n',...
%     (interfere_max- cons_sum_interfe_CF_FP_cog)>0);
% sum_rate_terres_CF_FP_cog = cal_sum_rate_terres(N,K,ww_CF_FP_cog,hGG,noise);
% fprintf('sum_rate_terres (alg_Closed_Form_FP) = %d (bps/Hz)\n',sum_rate_terres_CF_FP_cog);
% 
% 
% 
% %% Closed Form FP in cognitive network
% [ww_CF_FP_cog2,uu,vv,num_loop] = alg_Closed_Form_FP_cog2(N,K,MG,hGG,hGA,noise,pow_max,interfere_max);
% cons_sum_power_CF_FP_cog2 = cal_sum_power(N,K,ww_CF_FP_cog2);
% fprintf('%d/%d: the transmit power constraints of each terrestrial BS\n',...
%     sum((pow_max- cons_sum_power_CF_FP_cog2)>-1),length(pow_max));
% cons_sum_interfe_CF_FP_cog2 = cal_sum_interfe(N,K,ww_CF_FP_cog2,hGA);
% fprintf('%d/1: the interference power constraint of the primary network\n',...
%     (interfere_max- cons_sum_interfe_CF_FP_cog2)>0);
% sum_rate_terres_CF_FP_cog2 = cal_sum_rate_terres(N,K,ww_CF_FP_cog2,hGG,noise);
% fprintf('sum_rate_terres (alg_Closed_Form_FP) = %d (bps/Hz)\n',sum_rate_terres_CF_FP_cog2);















